function deflation_prob = posterior_predictive_distribution(WA,RA,QA,SA2,HZN,NMC);

% deflation_prob = posterior_predictive_distribution(WA,RA,QA,SA2,HZN,NMC);
% lnp is a sample from the posterior predictive density
% deflation_prob = probability that lnp < 0

% simulate log volatilities
Wchol = chol(WA,'lower');
e_vol_1 = Wchol(1,1)*randn(HZN,NMC);
e_vol_2 = (Wchol(2,1)/Wchol(1,1))*e_vol_1 + Wchol(2,2)*randn(HZN,NMC);
ln_r = log(RA(1)) + cumsum(e_vol_1);
ln_q = log(QA(1)) + cumsum(e_vol_2);
std_r = exp(ln_r).^.5;
std_q = exp(ln_q).^.5;
% simulate innovations to ln \mu, \pi
e_mu = std_q.*randn(HZN,NMC);
e_pi = std_r.*randn(HZN,NMC);
% future paths of \mu,\pi
mu = SA2 + cumsum(e_mu);
lnp = cumsum(mu + e_pi);
DD = zeros(size(lnp));
DD(lnp<0) = 1.0; 
deflation_prob = mean(DD,2);
end